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Quantum theory of the intra-band resonant tunneling of a Bose-Einstein condensate loaded in a 
two-dimensional optical lattice is considered. It is shown that the phenomena of quantum collapse 
and revival can be observed in the fully quantum problem. The meanfield limit of the theory is 
analyzed using the WKB approximation for discrete equations, establishing in this way a direct 
connection between the two approaches conventionally used in very different physical contexts. 
More specifically we show that there exist two different regimes of tunneling and study dependence 
of quantum collapse and revival on the number of condensed atoms. 



PACS numbers: 

I. INTRODUCTION 

Passage from the exact many-body description of an 
atomic gas at zero temperature to the mean-field the- 
ory is based on the assumption about large occupation 
number N of the ground (or, more generally, some spe- 
cific quantum) state. Then iV -1 / 2 can be shown to be 
the parameter of expansion, which in the leading order 
gives the Gross-Pitaevskii equation pQ . Respectively, the 
mentioned requirement must be also applied to mean- 
field theories of the Landau-Zener tunneling [2J [3J 0] , for 
which interesting experimental data are available [5], to 
instabilities of BEC's in optical lattices [3J E] observed 
experimentally in [7] , and to a theory of nonlinear Bloch - 
band tunneling [8 . Similar approach was also exploited 
in a pure quantum theory of a Bose-Einstein condensate 
(BEC) in a double-well trap [9], based on a model with 
linear coupling [10 . Either for the Bloch-band tunneling 
problem or for the tunneling in a double-well potential 
the approach of large N is justified when the system is 
close to the mean-field regime and the respective dynam- 
ical solutions are characterized by the essentially nonzero 
populations of the quantum states between which tunnel- 
ing occurs. Such regimes indeed exist when the descrip- 
tion is reduced to a dimer [TT|. At the same time the 
mean-field approximation to the above problems allows 
solutions where the populations of the states repeatedly 
become zero. Strictly speaking this violates the initially 
made supposition about large atomic numbers, because a 
state with negligible population cannot be treated in the 
mean-filed approximation. Therefore the theory requires 
modifications. 

We thus can formulate the main goal of the present pa- 
per as the analysis of the macroscopic nonlinear tunneling 
between two quantum states with the same (or infinitely 
close) energies in the limit corresponding to large total 
number of particles, allowing however populations of each 



of the states to become negligibly small at some moments 
of time. We will pay the main attention to the quantum 
effects, not accounted by the standard mean-field theory. 

To be specific, as a physical system we explore a BEC 
loaded in a two-dimensional optical lattice and consider 
tunneling between two states with the same energy. Such 
a situation can be experimentally realized in at least two 
different ways. First, one can use a non-separable lat- 
tice, with properly chosen parameters providing closer of 
the lowest gap (as, for example, this is suggested in [S]). 
Then nonlinear tunneling occurs between X and M points 
of two different bands (see Fig. [I]). More general situ- 
ation, however, corresponds to the nonlinear tunneling 
between two X points of the same band, rotated by n/2 
with respect to each other (these are the points Xi,2 in 
Fig. [I]). Such intra-band tunneling does not depend on 
the size of the gap, and hence can be observed in any 
lattice including separable one. Moreover, whenever one 
deals with the nonlinear tunneling, the stability issue ac- 
quires especial importance [3J [S] , since depending on the 
sign of the effective mass the homogeneous excitations 
can be either modulationally stable or unstable. This in 
particular means that tunneling between two states in 
different bands (like tunneling between X and M points, 
mentioned above, or, tunneling in one-dimensional case 
between two neighbor zones) represents a transition be- 
tween stable and unstable states, what causes asymmetry 
of the tunneling. That is why, in the present paper we 
concentrate on intra-band tunneling between two stable 
X points, what rules out instability and allows one to 
limit the consideration to plane Bloch waves. 

In the mean-field approximation this problem was con- 
sidered in |12j . which in the two-mode approximation is 
reduced to a model imitating cither a Josephson junc- 
tion or Rabi oscillations in a system of two-level atoms. 
Such a model describes population oscillations between 
two states. Using the mentioned analogy, one can predict 
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FIG. 1: (Color online) The contour plot of the periodic po- 
tential, given by equation ( |38[ ) below, used for the numeri- 
cal simulations. The inset shows the first Brillouin zone, the 
reciprocal lattice basis vectors Qi,2 and the high symmetry 
points F, M, and Xj. The nonlinear intra-band tunneling 
occurs between the points Xi and X2. 



that quantization of the motion should result in quan- 
tum collapse and revival of the above oscillations 1 3J . 
Obtaining those phenomena in the process of intra-band 
tunneling, constitutes the second aim of the present pa- 
per. 

Finally, we emphasize an interesting mathematical as- 
pect of the problem at hand. Namely, we will establish a 
formal link between the mean-field approximation, where 
the small parameter can be identified as 1/N and the 
discrete WKB method, having h as a formal small pa- 
rameter, which is described in details in Ref. [14] . In our 
approach the both methods, appearing originally as in- 
dependent as they are used in different setups and based 
on different physical assumptions, are intimately related 
to each other allowing us to appreciate the accuracy of 
the mean-field model, what will be the third goal of the 
present work. It is to be mentioned that earlier the link 
between the semi-classical and mean-field approaches was 
discussed in [TS], where a linearly couped two-level sys- 
tem, modeling either a BEC in a double-well trap or a 
spinor condensate, was considered and where the quan- 
tum corrections appeared as a decoherence of the quan- 
tum states. We, thus present one more physical system, 
a BEC in an optical lattice, having no linear coupling 
and rather different stability properties, which also al- 
lows experimental verification of the mean-field approx- 
imation. In our case the quantum effects will be mani- 
fested through the quantum collapse. 

The paper is organized as follows. In Sec. II we de- 
duce the quantum two-mode model describing resonant 
intra-band tunneling, expand its solution over the Fock 
basis and deduce the dynamical system for the expan- 
sion coefficients. In Sec. Ill we discuss the mean-field 
approximation starting with the derived dynamical sys- 
tem and based on the WKB method for discrete equa- 
tions. Sec. IV is devoted to numerical study of the both, 



full quantum and mean- field models. The outcomes are 
summarized in the Conclusion. 



II. QUANTUM MODEL 

A. Inter— and intra— band transitions 

Let us start with the Hamiltonian of a BEC in a two- 
dimensional, x £ M 2 , optical lattice V(x), 

H = J d 2 x^(x) (~^V 2 + V(x)^ 4>(x) 
v 

+| j d 2 x^ t (x)V t (x)V(x)V(x), (1) 

V 

where V = Mvq is the total area of the lattice consisting 
of M cells each one of the area vq, g is the interaction 
coefficient in two dimensions, and ip'(x) and f/>( x ) are 
the creation and annihilation field operators. Introduc- 
ing the Bloch waves <^„k( x ) through the standard eigen- 
value problem (-^V 2 + V(x) \ ip„ k = E n]i ip nk , where 

n is a number of the zone and k is the wave-vector in the 
first Brillouin zone (BZ), and considering them orthonor- 
malized, / v <£ nk (x)y> n / k /(x)dx = <W<W, we represent 



^( x ) = E Wk( x )b„ jk , 



(2) 



where the creation and annihilation operators satisfy 
the usual commutation relations [6„ k , 6 n 'k'] = and 

[b n k,b n , k ,] = <5rm'£>kk'- 

The expansion ([2| allows one to rewrite the Hamilto- 
nian in the form 

n,k 

Akik 2 k 3 k4 k l+ k 2-k3-k4,Q C 7i 1 ki ( ™2k 2 Un 3k3 Un 4k 4 i 

k 1 — k 4 

n 1: ...,n 4 

(3) 

where Q is an arbitrary vector of the reciprocal lattice, 
Xklk 2 2 k 3 3 k 4 4 = f / d2x <kX 2 k 2 ^3k3^uk 4 (4) 

V 

(hereafter an asterisk stands for complex conjugation). 

Next we apply the rotating wave approximation, i.e. 
neglect all four-wave processes which do not satisfy the 
energy conservation law. This allows us to simplify the 
Hamiltonian (J3j) : 

H = ^2 E nk bl k b nU + 

n , k 



+k 2 ,k 3 +k 4 +Q0,E„ lkl +B„ 2k2 ,E„3k 3 +E Tl4k4 



k!,...,k 4 

ni,...,n4 



y . v nin 2 n 3 n4i ) t (,t u u 

A A k x k 2 k 3 k 4 u n 1 ki u n 2 k 2 "™3 k 3 u n 4 k 4 



(5) 
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Further reduction of the Hamiltonian can be done by 
taking into account that the periodic medium is highly 
dispersive, and therefore for the resonant four-wave inter- 
actions, equality of the group velocities of the respective 
matter waves must be imposed. In the case of tunnel- 
ing between two states, say (ni,qi) and (71.2, q2), this 
means the constraint dE ni(il /dqi = dE n2 „ 2 /dq2- To- 
gether with the momentum and energy conservation laws, 
expressed by the first and the second Kronecker deltas in 
([5]) , this means that the resonant nonlinear tunneling can 
occur only between highly symmetric points of the BZ. 
As far as we have limited ourselves to the intra-band tun- 
neling, we conclude that the only points of the BZ the 
points X - satisfies the above requirements. Moreover, in 
the BZ there are only two physically different X-points. 



Assuming that the lattice is orthogonal with a period d 
in the both directions, the BZ can be identified simply 
as the domain [— 5 , f] X [— ^> f]- Consequently the two 
X-points we are interested in, below they are referred to 
as Xi and X2, correspond to the two vectors qi = (0, 5) 
andq 2 = (f,0). The group velocities in these points are 
zero. 

Now one readily verifies that ([5| is significantly simpli- 
fied allowing one to rewrite the Hamiltonian in the form 
of a sum 

H = J2 H n + H i-b (6) 

n 

where 




describes intra-band transitions between states Xi and X 2 of the n-th band and 



Hi-b 4 ^ ] Xqjqa (^niqi^n 2 q2 ^« 2 qi^niq2) (^niqi^«2q2 + ^n 2 qi^niq 2 ) ($) 

»1>2 



r 



describes transitions between X points of different bands. 
In the above formulas we introduced Xqjq 2 = Xkjk^kj 
and used the symmetry properties of the Bloch functions, 
giving in particular Xq? qi = Xq" q2 and Xq" q2 = Xq" qi > 
as well as the fact that y n qi and <p n q 2 correspond to X 
points and hence can be chosen real and periodic. 



B. Dynamical equations and their accuracy 

Dynamics of the condensate can be described by the 
time evolution of the coefficients of the expansion of a 
given multiparticlc state over the Fock basis. In 

the case when initially one band, say the band no, is pop- 
ulated more densely than the other ones, it is natural to 
represent the basis as \N%, N 2 ; n), where Nj stand for the 
occupation numbers of the states X^ of the given band, 
while n symbolically designate occupation numbers of all 
other bands. Thus 



No 



(9) 



N 1 ,N 2 ,n 



Now we consider the situation, where initially (i.e. at 
t — 0) all atoms belong to X points of the chosen band, 
i.e. when N\ + N 2 = N with N being the total number 
of particles. Then, introducing the simplified notation 
\Ni,N 2 ) = \Ni, N 2 ; 0) for such states, which also can be 
rewritten as \Ni,N 2 ) = \k, N — k) = \k)\N - k) (here 



k = N±), and respectively CN 1 N 2 o{t) 
express 



Ck(t), we can 



N 



N 



\m)=J2c k (0)\k,N-k), ^|C fc (0)| 2 = l. (10) 



fc=0 



fe=0 



For the next step, it is not difficult to verify that 
the dynamics of a state, initially spanned over the kets 
\k,N—k), i.e. subject to the initial condition (10 1, which 
is induced by the Hamiltonian ^ results only in the 
states belonging to initial sub-space, i.e. 



A'* 



N 



fc=0 



for all t. 

Indeed, substituting ^ in the Shrodinger equation 



.g|g) 



H\9) 



(12) 



and applying (iVi,iV2;n| in order to obtain the equa- 
tions for the expansion coefficients, one finds that for 
all coefficients with Ni + N 2 < N (i.e. for states with 
n^O) the derivatives dCN 1 N 2 n/dt linearly depend on dif- 
ferent Cjv< n 2 h' with N[ + N 2 < N, but do not depend on 
Cn 1 n 2 o- This follows from the relations Hi-b\Ni, N 2 ) = 
and H n \N 1} N2) = for n ^ n . The first of these 
formulas means that no coupling occurs when at least 
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one of the states in not occupied, while the second re- 
lation means zero result when one probes the energy of 
an empty band (H n does not originate inter-band transi- 
tions) and in our case the only band, no — 0, is initially 
occupied. Thus all the coefficients Cn 1 n 2 ii with n ^ 
are identically equal to zero, provided they are zero at 
t = 0. 

Summarizing the described situation subject to as- 
sumption that only one band no is populated and ne- 
glecting the unessential for the dynamics linear energy 
term -E„ oqi + E no<i2 , one arrives at the two-mode model 
whose Hamiltonian can be written down in the form 



H,, 



= Xu 



H 

2 



v noqi 



{n\ 



A 



(13) 



dimensional discrete space k = {0, 1, 2, . . . , N}. Indeed, 
setting x — k/N and h — 2/N, introducing a differen- 
tiable function ip(x) of the continuous variable x <G (0, 1), 
such that ip(x) = ^/NCk at the points x = k/N [5D] let 
us define the quantities: 



b<t\x) = b h (x + ft/2) ± b h (x - ft/2), 



b h {x) = [x(x + ft/2)(l - x)(l -x+ ft/2)] 1/2 , 
a(x) = 1 + 2(2A- l)bo(x). 



Now equation ( 14 I becomes 



where b 3 = 6„ oqj , i.e. b\\k,N-k) = y/k + l\k + l,N -k) 
and b\\k,N - k) = y/N - k + l\k,N - k + 1), rij = 
bn oClj bnoqj are the populations of the states Xj of the 
band n , Xu = X^"°, and A = Xq?q 2 ° /Xq?q? ■ Hamilto- 
nian ( 13 I preserves the total number of atoms in the Xj- 



points, what naturally reflects the approximations made. 

The above calculations, which resulted in ( 13 ), hold for 
the more general Hamiltonian (|5| obtained in the rotat- 
ing wave approximation, but fail for the original model 
accounting all possible transitions (not only the resonant 
ones). This defines the accuracy of the two- mode approx- 
imation: the ratios between the accounted and neglected 
terms are of the order of 1/N. 



Now the Schrodinger equation (12) with H — Ho re 



suits in a system of ordinary differential equations for the 
coefficients Cfc 



i dCfc 



A 



(frfc-lCfc-2 + bk+iCk+2) + yCl, 



(14) 



with the dimcnsionlcss time r = {Axi\N/K)t and with 
the coefficients 



a k = 1 + 2(2A - 1) 



1 - 



k fk 
N \ N 



1 

N 



N 

k 
N 



k 
N 



1 

N 



(15) 



(16) 



III. FICTITIOUS PARTICLE 
REPRESENTATION AND THE SEMICLASSICAL 
LIMIT 



ihd T ip{x) = ^-^-ip(x) 



ip(x). 



(17) 



Here the parameter ft plays the role of the Plank constant 
ft, p = ~ihd x is the "momentum operator" of the ficti- 
tious particle, and the "wave function" ^{x) is assumed 
to be normalized in the usual way L da:|^(a:)| 2 = 1. 



Equation ( 17 1 describes a fictitious quantum particle 



of the mass proportional to 1 /A moving in a compact 
curved space defined by the interval (0, 1) (hence b^\x) 
at the functions of the momentum operator) . Note that 
the "momentum" eigenvalues can be restricted to the first 
BZ, i.e. ' r 



d' di 



The semiclassical dynamics corresponds to the limit 
ft — > 0, i.e. when the number of BEC atoms N — > 00, 
what is the usual limit of the Gross-Pitaevskii equation. 
It should be noted that the characteristic time t of the 
evolution scales as (gN)~ 1 r, hence the quantity gN must 
be kept fixed, which is the second condition of the mean- 
field limit. It is also clear that the limit ft — > 0, if it 
exists, corresponds to the continuous limit of the discrete 
equation (17) (see also Ref. [14J). 



Mean-field equations viewed as a quasi-classical 
limit. 



Equation (14 1 can be viewed as the Schrodinger 



equation for a fictitious quantum particle in the one- 



In order to derive the quasi-classical equation corre- 
sponding to the limit ft — > we proceed in the usual 
way. Setting ip{x, r) = e lS (x,Tji)/h f or a conl pi ex action 
S(x, t, ft) viewed as a series S = + ftSW + wo 
get the equation 



Assuming that the function S(x, r, h) has derivatives 
with respect to x up to the second order, expanding 



S(x±h)-S(x) 



= ±4 0) + 0(h) 



and setting h — > in the equation ( 18 1 we get 
the Hamilton-Jacobi equation for the classical action 
S(x,t) = S(°\x,t)-t/2: 



S T = b, 



o(x){Acos(^) + 2A-l} 



(19) 



[recall that bo(x) = x(l — x)]. 

Thus the physical sense of the transition to the classi- 
cal limit in our problem is the transition from a discrete 
equation to its continuous limit. Note that in our setup 
the usual classical limit h — > 0, sometimes understood as 
a mathematical abstraction, acquires the well established 
sense of the limit of a large number of atoms and thus 
can be studied experimentally. 

In the case of a finite h (finite N) we have [p, x] = —ih, 
i.e. the usual canonical commutator of the momen- 
tum and coordinate. It turns out, however, that quasi- 
classical dynamics is more convenient to describe in terms 
of variables z — 1 — 2x and $ = S x = p, where x and 
p are the classical limits of the corresponding quantum 
variables. Then the Poisson brackets of the respective 
classical dynamics read 

{$, z} = lim % - [p, 1 - 2x] = -2. 

h^O h 

and the classical Hamiltonian can be recovered from the 
Hamilton-Jacobi equation ( fl9| ): 

1 



H = -(l- ^)(2A + Acos$ - 1). 



(20) 



Thus the classical Hamiltonian equations have the form 
z = -2^— = -(l-z 2 )sin$, (21a) 

$ = 2— = (1 -2A- Acos$)z. (21b) 

oz 

In our case these equations correspond to the mean field 
approximation . 

In order to clarify the physical meaning of the intro- 
duced classical variables z and <& let us compute 



(*|(n 2 -m)|*) = iV c k 



N 2 -N : 



N 



N 1 +N 2 = N 



N(l-2(x)), (22) 



and 



JV 



m{bl) 2 b\\^) = N 2 Y,Clb k+l C k+2 



k=0 



= N\b h [x+-)e^) (23) 



(in the last formula we used that Ck+2 = e* p Cfc). 
Thus in the limit h — > 



z = lim (z) 



(n 2 ) - (rti) 
N 



(24) 



is the relative population of the states (as it is clear z S 
[-1,1]) and 



*= lim(p)=arg(^) 2 6?). 

h—>0 



(25) 



is the relative phase. 



B. On the mean-field dynamics 

Let us now discuss some aspects of the dynamics de- 
scribed by the mean-field approximation. First of all 



we observe that the system (21 1 has two fixed points: 



Pi 



0,$ 



0) and P 2 = (z = 0,$ = tt). 



The respective frequencies are Q\ = 4^ (A — A c ) and 
fl 2 — ^(1 — A) with A c = 1/3. Thus one easily verifies 
that P 2 is a local maximum, while Pi is a saddle-point 
for A < A c and a local minimum for A > A c . 

Next, following Ref. [TJ] we take into account that the 
classical energy E — H(z, $) is bounded by the two po- 
tential curves, one corresponding to <i> = [the curve 
J7 (+ )(a;)] and the other to <I> = n [the curve [A - ) (a;)]: 

£/(-)(*) = (A - l)b (x), U(+\x) = (3A - 1)60(1), (26) 

i.e. U(~\x) < E < U^ + \x). Therefore, the turning 
points of the classical system ( |21a l-( 21b I lie on the curves 
U^\x). One can express the period of oscillations be- 
tween two turning points X\ and x 2 as follows 



T 



Ax 



y/(UM -E){E -[/■(->)" 



(27) 



This integral can be expressed in terms of the complete 
elliptic integral of the first kind K(-). To this end we 
single out four physically different cases: 
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Case 1: If A < A c and the energy of "classical" motion 
is (3A — l)/4 < E < 0, the turning points belong to 
different curves, for instance: 



J'2 



1 - \ 1 



AE 
1 - A 



1-1/1 



AE 
1 -3A 



[/(+) 



and the period is given by 



AK 



2y/-2gA 



- V /(1-3A)(1-A+4B) 

v / (l-3A)(l-A + 4£) 



(28) 



The difference in the atomic population oscillates about 
a nonzero value. 

Case 2: If A < A c and the energy of "classical" motion 
is (A — I) /A < E < (3A - 1)/ A < 0, the turning points 
belong to a single curve: 



1 



Xl,2 



l±\ 1 



4£ 
1 - A 



and the period is given by 
2V2 



T 2 = K 



V(l-3A)(l-A + 4g) 
2V-2.EA 



(29) 



The difference in the atomic population oscillates about 
zero. 

Case 3: If A > A c and the energy of "classical" motion 
is in the interval (A — l)/4 < E < 0, the turning points 
belong to a single curve: 




V /(1-A)(3A-1-4S) 



(30) 



The two states have equal average populations. 

Case 4- If A > A c and the energy of "classical" motion 
is in the interval < E < (3A — 1) /A, the turning points 
belong to a single curve: 



1 



Xl,2 



1±W1 



AE 



2 V V 1 - 3A 
and the period is given by 



8K 



1+^ 



v/(3A-l)(l-A-4£) 



(31) 



The two states have equal average populations. 

Finally we mention that equations (21 1 can be directly 



obtained from the two-mode Hamiltonian of the Bloch- 
band tunneling in the meanfield approximation (see |12l 
[TH] for more details): 



H, 



m—f 



Xn\Ai 



+ Xii\A 2 \ 4 + Axi2\A 1 \ 2 \A 2 \ 2 
-xMA^Al + xMA^Aj. (32) 



where the complex amplitudes A1.2 are determined by 
the expansion of the order parameter ^ 



# = ./4i(f)<Pi(r)e- 



iEI 



A 2 (t)<p 2 (r)e 



-iEt 



(33) 



and determine average populations of the levels JVj = 
\Aj\ 2 , and <Pj(r) are the respective Bloch states. The 
complex amplitudes are connected by the particle conser- 
vation law: N = N\ + M 2 . Now the system (pTal , ((21b I 



is obtained from the Hamiltonian equations by defining 
z JA^_W^_ md $ = argM?( ^ )2} (34) 



[c.f. Eq. (251] 



C. "Coherent" states 

Turning now to the quantum system, we observe that 
the dynamical equations for the coefficients C n have to 
be supplied by the initial conditions. As in the stan- 
dard WKB approximation, the corresponding initial con- 
ditions must be smooth enough. At the same time a natu- 
ral question arises about construction of quantum states, 
\^(t)) c , most closely resembling the meanfield dynamics. 
We will refer to such states as coherent states. 

In order to construct such states we recall the ex- 



plicit form of the meanfield ansatz ( 33 1 and consider 
the respective boson operator 



aib[ + a 2 b 2 , where 
ax 2 are time dependent complex parameters satisfying 



a 2 



1. Next we define 
(ai(i)fcj 



\oti,a 2 ) 



a 2 {t)b\) N 



0,0) 




N\ 



k\(N-k)\ 



al<x%- k \k,N -ft). (35) 



By using the identities: 



(a 1 ,a2\n j \ai,a 2 ) = ^l^l 2 , 

(a l ,a2\n 2 \a 1 ,a 2 }=N\a j \ 2 +N(N-l)\a j \ 4 , 

(ai, a 2 \nin 2 \ai, a 2 ) = N(N - l)|aia 2 | 2 , 

(a 1 ,a2\{bl) 2 b 2 1 \a 1 ,a 2 }=N(N-l)a 2 1 {a* 2 ) 2 , 

we obtain for the energy (dropping an inessential con- 
stant) 

E={a\H\a)=xuN(N-l) 

x {K| 4 + M 4 + A[4M 2 |a 2 | 2 + a\{a* 2 ) 2 + a 2 (al) 2 ]} . 

(36) 
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Writing the Hamilton equations for the complex ampli- 
tudes a-j in the form 



dt 



1 dE 



in terms of the normalized time r = A\nN /h and setting 



"2 



|«i| 2 and $=axg{a\af} (37) 



[c.f. Eq. (34 1] we recover the system (211 in the limit 
N — > oo. Comparing (32 1 with (36) one gets the physical 



meaning and a link among the quantum mechanical and 
mean-field amplitudes <x, and Aj. 



IV. NUMERICAL SIMULATIONS 

In order to proceed with the numerical simulation of 
the described phenomena we observe that their occur- 
rence does not depend on a particular choice of the po- 
tential. Therefore we concentrate on the simplest case of 
a lattice having cos-like profile. Recalling that the spa- 
tial coordinates are measured in the units of the lattice 
period d, the energy is measured in terms of the recoil 
energy E r = h 2 -K 2 /(2mcP), the time is measured in the 
units H/E r and the lattice is 7r-periodic, we set 



V = Vb[cos(2a:) + cos(2y)]. 



(38) 



The respective BZ is given by [—1,1] x [—1,1]. 

We consider intra-band tunneling between two X- 
points of the second lowest band. The respective tensors 
of the inverse effective mass have positive components, 
and therefore the respective Bloch states are modula- 
tionally stable, provided the interactions among atoms 
are repulsive (i.e. Xij > 0). 

Generally speaking one may have two physically dif- 
ferent situations. When Vo > Vth, with Vth ~ 0.627 for 



the case ( 38 ) , there exists a full two-dimensional gap in 



the spectrum, whose width depends on a particular value 
of Vq- If the potential amplitude is weak enough, i.e. 
Vo < Vth j the gap is closed (the gap width becomes zero 
at Vq = Vth). For Vo = Vth one computes, using (38 1, 
that m 0.1494 and x q 2 q2 ~ 0.1303, and respec- 

tively A = Ath ~ 0.8727 > A c . Thus, generally speaking, 
one can distinguish three different regions of the param- 
eters, which correspond to (i) A < A c , (ii) A c < A < A t h, 
and (iii) A > A t h- The parameters A t h and A c have how- 
ever different physical origins: Ath is related to the band 
structure and thus can affect physical applicability of the 
two mode model, due to possible tunneling to the other 
bands, while A c is an intrinsic parameter of the model. 

Since our numerical study aims to check precisely the 
two-mode model we will consider only the case where the 
full gap is open, i.e. V > V t h, and select two particular 
cases: the Case I with Vb = 6 when X q ^ qi ~ 0.2708 and 
X q 2 iq2 « 0.0.0652, and respectively A « 0.2407 < A c , 



and the Case II with Vq = 1 when Xqiqi ~ 0.1551 and 
X q 2 q2 « 0.0626, and respectively A « 0.4040 > A c . 

The numerical simulations of the discrete Schrodinger 
equation ( 14 1 were performed by using the variable order 
Adams-Bashforth-Moulton solver. The error was con- 
trolled by checking the norm of the vector Ck(t). The 
values of z and $ were obtained by using the correspon- 
dence formulas ( 22 1 and ( 23 1 . 



In Fig. [2] we show typical dynamics of the populations 
in the Case I, i.e. for A < A c , of the condensate of N — 
350 atoms. One can observe relatively fast, i.e. at t ~ 
100, decay of Rabi oscillations followed by almost steady 
distribution of the atoms, approximately 90% of atoms 
concentrated in one state. This is behavior is typical 
for a quantum collapse. The both populated states are 
characterized by the same phases. After much longer 
interval of time At ~ 1500 the revival of the oscillations 
is observed. 
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FIG. 2: Quantum evolution for N = 350 BEC atoms and the 
lattice with A = 0.2407. As the initial population we used 
the Gaussian distribution Cfc(0) = exp{ifc — (fc — fe c ) 2 /(2cr 2 )} 
with a 2 — 50 and k c — 50. 



Although collapse of Rabi oscillations was earlier ob- 
served in pure mean-field models (sec [3, 8 ), the na- 
ture of the phenomenon considered here is very different. 
Suppression of the oscillations in the case of spatially 
extended systems is related to development of inhomo- 
geneous spatial patterns, mainly related to modulational 
instability of one of the states. In the case at hand, how- 
ever the effect is essentially quantum and disappears in 
the meanfield limit (where N — > 00). This is clearly il- 
lustrated in Fig. [3] where we compare the results of the 
quantum dynamics with its mean-field limit at earlier 
stages of the evolution. 

Now we turn to the Case II, where A > A c , which is 
shown in Figs. [4] and [5] The main feature of this situation 
is that as a result of quantum collapse both states become 
equally populated (z — > 0), what is also in accordance 
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FIG. 3: The initial quantum evolution (solid lines) vs mean- 
field evolution (dashed lines) for parameters as in Fig. [5] 



with the meanfield dynamics (c.f. Figs. [2] and [3]). The 

revival occurs at latter times At ~ 3000. 
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FIG. 4: Quantum evolution for N = 500 BEC atoms and 
the lattice with A = 0.404. The initial condition used is 
Cfc(0) = exp{-(fc - k c ) 2 /(2a 2 )} with k c = 25 and a 2 = 50. 
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FIG. 5: Quantum evolution for N = 1000 BEC atoms and 
the lattice with A = 0.404. The initial condition used is 
C fc (0) = exp{ifc - (jfe - k c ) 2 /{2a 2 )} ($(0) = 2) with k c = 10 
and a 2 = 50. 



performed at a constant gN, i.e. at a constant effective 
nonlincarity of the system. Meantime, the tunneling time 
depends on the relative value of the nonlinear interaction 
term, i.e. on A. 

Comparing Figs. [2] and [4] one observes the expected 
effect of the delay of the quantum collapse in terms of 
the dimensionless time r with increase of the number of 
particles. At the same time the numerical simulations 
did not reveal any significant effect of the initial phase 
mismatch $(0) on the collective dynamics. 

Another relevant parameter is the initial distribution 
of the atoms. Clearly, the closer is the initial distribution 
of atoms (defined by k c and a in the numerics) to the one 
violating the mean- filed assumptions, when almost all of 
the atoms are at one of the Xj points, the further is the 
dynamics from the quasi-classical one, which is reflected 
in Fig. [5] in the fact that the collapse time of the oscil- 
lations is twice as smaller as that of Fig. [4] for twice as 
much number of atoms. Also, in Fig. [5] the phase dynam- 
ics at the recurrence of quasi-classical relative population 
oscillations contains the quasi-classical oscillations inter- 
rupted by the 2irm jumps. 



In all the figures one can observe that the period of 
the Rabi oscillations is very accurately reproduced by 



the mean-filed approximation, i.e. by the formulas (28 1, 
|29j), (f30j), and (|30|> . In particular, Figs. [| (or^2| 
and [5 correspond to the Cases 1 ( with z = 0.72 and 
E fa -0.033) and 3 (with z = 0.98 and E ps -0.0059) 



described in Sec. IIIB Then using (28 1 and (3l| one 
computes T pa 16.86 and T ps 51.99, what matches very 
well with the periods obtained from the direct numerical 
simulations. The described behavior has simple physi- 
cal explanation: passage to the mean-field description is 



V. DISCUSSION AND CONCLUSION 

In the present paper we have considered the quantum 
tunneling of a Bose-Einstein condensate loaded in a two 
dimensional lattice. The considered tunneling occurs be- 
tween two stable states, hence allows us to restrict the 
consideration to a spatially homogeneous model. Such 
a model, in the rotating wave approximation, was fur- 
ther reduced to the effective two-state system, where 
both states possess the same energies. The considered 
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tunneling is essentially nonlinear and related to the four 
wave-mixing due to the two-body interactions (in a pure 
linear system such a tunneling cannot occur because of 
violation of the momentum conservation). 

The main finding of the work is the quantum collapse 
and revivals of the Rabi oscillations between the tunnel- 
ing states, which are suppressed in the mean-field ap- 
proximation due to the fact that tunneling between two 
equal-energy states result in configurations of low atomic 
population of one of the states and thus, strictly speak- 
ing, cannot be described in terms of the mean-field ap- 
proximation. Nevertheless, the latter approach turns out 
to be useful in predicting different regimes of the tunnel- 
ing, with collapse occurring to either equal or disbalanced 
populations of the states, and for accurate estimate of 
the frequency of the Rabi oscillations (i.e. the tunneling 
time). 

Existence of a well-defined frequency of the oscillations 
and the difference in spatial patterns of the condensate 
in different X-states of the lattice suggest a way of exper- 
imental observation of the phenomenon, which could be 
based, for example, on direct imaging at specific moments 
of time. To estimate the physical time scale correspond- 
ing to the period of the Rabi oscillations, we consider 
a condensate of 87 Rb atoms, with the s-wave scattering 
length a s = 5.25nm, which is tightly trapped, say in 
y direction, with the respective linear oscillator length 
£ = 0.1 /im, and has N — 350 atoms occupying M = 25 
sites of a lattice with the period d = 2 /im (i.e. the 2D 
optical lattice is imposed in the (x, z) plane and has 5 
cells in each direction). In this case g (the 2D nonlinear 
constant) is given by g = 2y/2irh 2 a s /m£ y . Then using 
the data from Fig. [2] (also Fig. |3| we obtain that the 
dimensionless period T ~ 16.68 in the physical units is 
Tgb ~ 0.00655 s. Thus, taking into account that the col- 
lapse and subsequent revival would occur after about 1 s 
and 10 s, respectively, we see that this is a relatively slow 
process. This is an expected slowness as far as in this 
example we used a small number of atoms not providing 



effective enough tunneling due to the two-body interac- 
tions. The collapse and revival time could be reduced by 
using a larger condensate, however, for a very large num- 
ber the collapses and revivals will be suppressed because 
of approaching the mean-field limit. An alternative way 
to accelerate the quantum collapse, and thus make easier 
its observation is to use more light atoms, say sodium 
ones, or by increasing the scattering length. 

A number of important questions, however, are left for 
further studies. Among them we mention the resonant in- 
teraction of more than two states, the quantum theory of 
nonlinear Landau-Zener tunneling, the study of interplay 
of quantum collapse and of modulational (dynamical) in- 
stability in the case when the spatial extension of the 
system is taken into account and the states among which 
the tunneling occurs possess different stability properties 
(what is an intrinsic feature of the nonlinear systems), 
the effect of opening or closer of the total gap in the 
lattice, the interplay between the quantum collapse and 
the dynamical collapse (also called the blow up) which 
can occur in Bose-Einstein condensates with attractive 
atomic interactions, etc. 
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